function [Q,grad,se,p1mat] = func_mixed_logit(theta,y,X,draws,sumidx,seoption)

beta = theta(1:size(X,2));
sigma = theta(size(X,2)+1:end);


sumdraws = zeros(size(draws,1),size(draws,3));
for k=1:size(draws,2)
   sumdraws = sumdraws + sigma(k)*squeeze(draws(:,k,:));
end

umat = (X*beta)*ones(1,size(draws,3)) + sumdraws;

p1mat = exp(umat)./(1+exp(umat));
p0mat = 1-p1mat;

ymat = y*ones(1,size(draws,3));

probmat = p1mat.*ymat+p0mat.*(1-ymat);

Lmat = zeros(max(sumidx),size(draws,3));

for s=1:size(draws,3)
    Lmat(:,s) = accumarray(sumidx,probmat(:,s),[],@prod);
end

countidx = accumarray(sumidx,ones(size(sumidx)),[],@sum);

Lmatmean = mean(Lmat,2);

Q = -mean(log(Lmatmean(countidx>0)));

inner = (2*ymat-1).*(1-probmat).*Lmat(sumidx,:);

grad = zeros(size(theta));

for b=1:length(beta)
    tempgrad = accumarray(sumidx,X(:,b).*mean(inner,2),[],@sum)./Lmatmean;
    grad(b) = -mean(tempgrad(countidx>0));
end

for b=length(beta)+1:length(theta)
    tempgrad = accumarray(sumidx,mean(squeeze(draws(:,b-length(beta),:)).*inner,2),[],@sum)./Lmatmean;
    grad(b) = -mean(tempgrad(countidx>0));
end

if seoption==0
se = [];
else
varcov = zeros(size(theta,1));  

fderiv = zeros(max(sumidx),size(theta,1));

for b=1:length(beta)
    fderiv(:,b) = -accumarray(sumidx,X(:,b).*mean(inner,2),[],@sum)./Lmatmean;
end

for b=length(beta)+1:length(theta)
    fderiv(:,b) = -accumarray(sumidx,mean(squeeze(draws(:,b-length(beta),:)).*inner,2),[],@sum)./Lmatmean;
end

for b=1:size(theta,1)
    for c=1:size(theta,1)
        varcov(b,c) = mean(fderiv(countidx>0,b).*fderiv(countidx>0,c));
    end
end

varcov = inv(varcov)/sum(countidx>0);



se = sqrt(diag(varcov));
end


